%Script to validate teh candidate models for IGRF12

%% plot data and secular variation

S = dir('/Users/manojnair/data/obs_mag_data/BGS data/*.nc');


% selected_obs = ['AAA';	'ABG';	'ABK';	'ASP';	'BDV';	'BEL';	'BFO';	...
%     'BLC';	'BMT';	'BOX';	'BOU'; 'CBB';	'CLF';	'CTA';	'CZT';	'DRV';	...
%     'EBR';	'EYR';	'FCC';	'FUR';	'GNA'; 'GLN';	'HLP';	'HRN';	'KAK';	...
%     'KDU';	'KNY';	'LRM';	'LVV';	'MAW';	'MBO';	'MCQ';	'MEA';	...
%     'MMB';	'NGK';	'NUR';	'NVS';	'OTT';	'PAF';	'RES';	'SOD';	...
%     'STP';	'STJ';	'THY'; 'TUC';	'UPS';	'VIC'];

selected_obs = cell2mat(textread('/Users/manojnair/projects/wmm2015_validation/selected_obs.txt','%s'));

time_array_hour = (datenum(2010,1,1,0,0,30): (1/(24)): datenum(2014,12,31,23,59,30))';

for i = 1:length(S),
    filestring(i,:) = S(i).name(1:3);
end;

[data_obs, ia, ib ] = unique(filestring,'rows');

[spline_obs, ia, ib]= intersect(data_obs, selected_obs, 'rows');

% spline_obs = filestring; % using all observatories

% plot

f = figure(1);
set(f,'Position',[40          19        1557         937]);


annual_time_1 = datenum(2010,1:48,15);
annual_time_2 = datenum(2011,1:48,15);
actual_sec_time = datenum(2010,7:54,17);
annual_time = datenum(2010,1:60,15);

igrf_A = '/Users/manojnair/projects/igrf12_validation/SV-2015-2020-A.txt';
igrf_B = '/Users/manojnair/projects/igrf12_validation/SV-2015-2020-B.txt';
igrf_C = '/Users/manojnair/projects/igrf12_validation/SV-2015-2020-C.txt';
igrf_D = '/Users/manojnair/projects/igrf12_validation/SV-2015-2020-D.txt';
igrf_E = '/Users/manojnair/projects/igrf12_validation/SV-2015-2020-E.txt';
igrf_F = '/Users/manojnair/projects/igrf12_validation/SV-2015-2020-F.txt';
igrf_GB = '/Users/manojnair/projects/igrf12_validation/SV-2015-2020-GB.txt';
igrf_HB = '/Users/manojnair/projects/igrf12_validation/SV-2015-2020-HB.txt';
%igrf_I = '/Users/manojnair/projects/igrf12_validation/SV-2015-2020-I.txt';
% replacing the I model with IGRF11 for a quick check.
igrf_I = '/Users/manojnair/projects/igrf12_validation/igrf11.txt';


for i = 1:length(spline_obs),
    
    ncfname = ['/Users/manojnair/data/obs_mag_data/BGS data/' spline_obs(i,:) '_1995_2014.nc'];
    
    
    %display(ncfname);
    
    [x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj] = read_geomag_netcdf(ncfname, 131496, 43824, 0);
    
    if any(x_data),
        
        
        
        % fit splines
        L = isnan(x_data);
        
        data_array = x_data(~L);
        first_index = find(~L, 1, 'first');
        
        last_index = find(~L, 1, 'last');
        
        time_array = time_array_hour(~L);
        
        b1 = min(time_array):365:max(time_array)+10;
        
        if length(b1) > 1,
            
            sp_x=spline(b1,x_data(~L)'/spline(b1,eye(length(b1)),time_array'));
            sp_y=spline(b1,y_data(~L)'/spline(b1,eye(length(b1)),time_array'));
            sp_z=spline(b1,z_data(~L)'/spline(b1,eye(length(b1)),time_array'));
            
            if  obj.geospatial_lon > 180,
                obj.geospatial_lon = obj.geospatial_lon - 360;
            end;
            
            
            
            geoc_lat = geod2geoc(obj.geospatial_lat, 0, 'WGS84'); % geocentric lat
            geoc_r = ( geocradius(geoc_lat,'WGS84') - 6371.2*1000 )/1000; % geocentric alt in km
            
            % secular variation is linear and same through out the model so
            % no need to calculate it through the loop.
            
            
            geoc_lat = geod2geoc(obj.geospatial_lat, 0, 'WGS84'); % geocentric lat
            geoc_r = ( geocradius(geoc_lat,'WGS84') - 6371.2*1000 )/1000; % geocentric alt in km
            
            obs_data(i).name = spline_obs(i,:);
            
            obs_data(i).sec_A = magsynth(geoc_r,geoc_lat, ...
                obj.geospatial_lon,2015.0,igrf_A);
            obs_data(i).sec_B = magsynth(geoc_r,geoc_lat, ...
                obj.geospatial_lon,2015.0,igrf_B);
            obs_data(i).sec_C = magsynth(geoc_r,geoc_lat, ...
                obj.geospatial_lon,2015.0,igrf_C);
            obs_data(i).sec_D = magsynth(geoc_r,geoc_lat, ...
                obj.geospatial_lon,2015.0,igrf_D);
            obs_data(i).sec_E = magsynth(geoc_r,geoc_lat, ...
                obj.geospatial_lon,2015.0,igrf_E);
            obs_data(i).sec_F = magsynth(geoc_r,geoc_lat, ...
                obj.geospatial_lon,2015.0,igrf_F);
            obs_data(i).sec_GB = magsynth(geoc_r,geoc_lat, ...
                obj.geospatial_lon,2015.0,igrf_GB);
            obs_data(i).sec_HB = magsynth(geoc_r,geoc_lat, ...
                obj.geospatial_lon,2015.0,igrf_HB);
            obs_data(i).sec_I = magsynth(geoc_r,geoc_lat, ...
                obj.geospatial_lon,2015.0,igrf_I);
            sec_x_obs =  ppval(sp_x,annual_time_2) - ppval(sp_x,annual_time_1) ;
            sec_y_obs =  ppval(sp_y,annual_time_2) - ppval(sp_y,annual_time_1) ;
            sec_z_obs =  ppval(sp_z,annual_time_2) - ppval(sp_z,annual_time_1);
            %
            
            obs_data(i).sec_x_obs =  ppval(sp_x,datenum(2014,10,19)) - ppval(sp_x,datenum(2013,10,19)) ;
            obs_data(i).sec_y_obs =  ppval(sp_y,datenum(2014,10,19)) - ppval(sp_y,datenum(2013,10,19)) ;
            obs_data(i).sec_z_obs =  ppval(sp_z,datenum(2014,10,19)) - ppval(sp_z,datenum(2013,10,19));
            
            fprintf('%s \n', spline_obs(i,:));
            
            
            obs_data(i).lat = obj.geospatial_lat;
            obs_data(i).lon = obj.geospatial_lon;
            
            
            %             % remove non data region.
            %
            %             sec_x_obs(actual_sec_time > time_array_hour(last_index)) = NaN;
            %             sec_y_obs(actual_sec_time > time_array_hour(last_index)) = NaN;
            %             sec_z_obs(actual_sec_time > time_array_hour(last_index)) = NaN;
            %
            %
            %
            %
            %             f = figure(1);
            %             set(f,'Position',[ 560    39   778   909]);
            %             subplot(321);
            %
            %             plot(actual_sec_time, sec_x_obs,'LineWidth',2);
            %             set(gca,'FontSize',16);
            %             hold on;
            %             plot(datenum(2015,1,1), obs_data(i).sec_A(1),'r*','MarkerSIze',20);
            %             plot(datenum(2015,1,1), obs_data(i).sec_B(1),'r*','MarkerSIze',20);
            %             plot(datenum(2015,1,1), obs_data(i).sec_C(1),'r*','MarkerSIze',20);
            %             plot(datenum(2015,1,1), obs_data(i).sec_D(1),'r*','MarkerSIze',20);
            %             plot(datenum(2015,1,1), obs_data(i).sec_E(1),'r*','MarkerSIze',20);
            %             plot(datenum(2015,1,1), obs_data(i).sec_F(1),'r*','MarkerSIze',20);
            %             plot(datenum(2015,1,1), obs_data(i).sec_GB(1),'r*','MarkerSIze',20);
            %             plot(datenum(2015,1,1), obs_data(i).sec_HB(1),'r*','MarkerSIze',20);
            %             plot(datenum(2015,1,1), obs_data(i).sec_I(1),'r*','MarkerSIze',20);
            %
            %             % plot(datenum(2014,7,2), ppval(sp_x,datenum(2014,7,2)),'g*','MarkerSIze',20);
            %
            %
            %             h=legend('OBS','WMM2010','WMM2015NGDC','WMM2015BGS');
            %             set(h,'FontSize',8,'Location','NorthWest');
            %
            %
            %             datetick('x','keeplimits');
            %             ylabel('X(nT/yr)');
            %             title(sprintf('%s lat %5.2f lon %5.2f ',spline_obs(i,:), ...
            %                 obj.geospatial_lat, obj.geospatial_lon));
            %
            %             subplot(323);
            %             plot(actual_sec_time, sec_y_obs,'LineWidth',2);
            %             set(gca,'FontSize',16);
            %             hold on;
            %             plot(datenum(2010,1,1), sec_wmm2010(2),'r*','MarkerSIze',20);
            %             plot(datenum(2015,1,1), sec_wmm_ngdc(2),'k*','MarkerSIze',20);
            %             plot(datenum(2015,1,1), sec_wmm_bgs(2),'b*','MarkerSIze',20);
            %             %plot(datenum(2014,7,2), ppval(sp_y,datenum(2014,7,2)),'g*','MarkerSIze',20);
            %
            %
            %             datetick('x','keeplimits');
            %             ylabel('Y(nT/yr)');
            %
            %             subplot(325);
            %             plot(actual_sec_time, sec_z_obs,'LineWidth',2);
            %             set(gca,'FontSize',16);
            %             hold on;
            %             plot(datenum(2010,1,1), sec_wmm2010(3),'r*','MarkerSIze',20);
            %             plot(datenum(2015,1,1), sec_wmm_ngdc(3),'k*','MarkerSIze',20);
            %             plot(datenum(2015,1,1), sec_wmm_bgs(3),'b*','MarkerSIze',20);
            %             % plot(datenum(2014,7,2), ppval(sp_z,datenum(2014,7,2)),'g*','MarkerSIze',20);
            %
            %
            %             datetick('x','keeplimits');
            %             ylabel('Z(nT/yr)');
            %
            %
            %
            %             subplot(322);
            %
            %             plot(time_array_hour,x_data);
            %             hold on;
            %             spline_fit = ppval(time_array_hour,sp_x);
            %             spline_fit(1:first_index) = NaN;
            %             spline_fit(last_index:end) = NaN;
            %             plot(time_array_hour,spline_fit,'r');
            %             plot(datenum(2010:0.1:2014.9,0,0), wmm_2015_ngdc(:,1),'r.-','MarkerSize',15);
            %             plot(datenum(2010:0.1:2014.9,0,0), wmm_2015_bgs(:,1),'b.-','MarkerSize',15);
            %
            %             h=legend('OBS','FIT','WMM2015NGDC','WMM2015BGS');
            %             set(h,'FontSize',8);
            %             set(gca,'FontSize',16);
            %
            %             ylabel('X nT');
            %             axis([datenum(2010,1,1), datenum(2015,2,1), -inf inf]);
            %             datetick('x','keeplimits');
            %
            %
            %             subplot(324);
            %
            %             plot(time_array_hour,y_data);
            %             hold on;
            %             spline_fit = ppval(time_array_hour,sp_y);
            %             spline_fit(1:first_index) = NaN;
            %             spline_fit(last_index:end) = NaN;
            %             plot(time_array_hour,spline_fit,'r');
            %
            %             hold on;
            %             plot(datenum(2010:0.1:2014.9,0,0), wmm_2015_ngdc(:,2),'r.-','MarkerSize',15);
            %             plot(datenum(2010:0.1:2014.9,0,0), wmm_2015_bgs(:,2),'b.-','MarkerSize',15);
            %
            %
            %             set(gca,'FontSize',16);
            %             ylabel('Y nT');
            %             axis([datenum(2010,1,1), datenum(2015,1,1), -inf inf]);
            %             datetick('x','keeplimits');
            %             subplot(326);
            %
            %             plot(time_array_hour,z_data);
            %             hold on;
            %             spline_fit = ppval(time_array_hour,sp_z);
            %             spline_fit(1:first_index) = NaN;
            %             spline_fit(last_index:end) = NaN;
            %             plot(time_array_hour,spline_fit,'r');
            %
            %             hold on;
            %             plot(datenum(2010:0.1:2014.9,0,0), wmm_2015_ngdc(:,3),'r.-','MarkerSize',15);
            %             plot(datenum(2010:0.1:2014.9,0,0), wmm_2015_bgs(:,3),'b.-','MarkerSize',15);
            %
            %
            %             set(gca,'FontSize',16);
            %             ylabel('Z nT');
            %             axis([datenum(2010,1,1), datenum(2015,1,1), -inf inf]);
            %             datetick('x','keeplimits');
            %
            %
            %             set(f,'PaperPositionMode','auto');
            %             saveas(f,['/Users/manojnair/projects/wmm2015_validation/' sprintf('%s',spline_obs(i,:)) '_Sec_Var_NGDC_vs_BGS'],'epsc');
            
            
            %pause;
            
            clear  sp_x sp_y sp_z;
            %clf;
            
        end;
        
    end;
    
end;


%% extract data

for i = 1:length(obs_data),
    sec_A_x(i) = obs_data(i).sec_A(1);
    sec_A_y(i) = obs_data(i).sec_A(2);
    sec_A_z(i) = obs_data(i).sec_A(3);
end;

for i = 1:length(obs_data),
    sec_B_x(i) = obs_data(i).sec_B(1);
    sec_B_y(i) = obs_data(i).sec_B(2);
    sec_B_z(i) = obs_data(i).sec_B(3);
end;

for i = 1:length(obs_data),
    sec_C_x(i) = obs_data(i).sec_C(1);
    sec_C_y(i) = obs_data(i).sec_C(2);
    sec_C_z(i) = obs_data(i).sec_C(3);
end;

for i = 1:length(obs_data),
    sec_D_x(i) = obs_data(i).sec_D(1);
    sec_D_y(i) = obs_data(i).sec_D(2);
    sec_D_z(i) = obs_data(i).sec_D(3);
end;

for i = 1:length(obs_data),
    sec_E_x(i) = obs_data(i).sec_E(1);
    sec_E_y(i) = obs_data(i).sec_E(2);
    sec_E_z(i) = obs_data(i).sec_E(3);
end;

for i = 1:length(obs_data),
    sec_F_x(i) = obs_data(i).sec_F(1);
    sec_F_y(i) = obs_data(i).sec_F(2);
    sec_F_z(i) = obs_data(i).sec_F(3);
end;

for i = 1:length(obs_data),
    sec_GB_x(i) = obs_data(i).sec_GB(1);
    sec_GB_y(i) = obs_data(i).sec_GB(2);
    sec_GB_z(i) = obs_data(i).sec_GB(3);
end;

for i = 1:length(obs_data),
    sec_HB_x(i) = obs_data(i).sec_HB(1);
    sec_HB_y(i) = obs_data(i).sec_HB(2);
    sec_HB_z(i) = obs_data(i).sec_HB(3);
end;

for i = 1:length(obs_data),
    sec_I_x(i) = obs_data(i).sec_I(1);
    sec_I_y(i) = obs_data(i).sec_I(2);
    sec_I_z(i) = obs_data(i).sec_I(3);
end;


%% print ms results

sv_models = {'A','B','C','D','E','F','GB','HB','I'};
for i = 1:9,
    
    eval(['sec_x = sec_' cell2mat(sv_models(i)) '_x;']);
    eval(['sec_y = sec_' cell2mat(sv_models(i)) '_y;']);
    eval(['sec_z = sec_' cell2mat(sv_models(i)) '_z;']);
    
fprintf('%s %5.2f %5.2f %5.2f\n', cell2mat(sv_models(i)),rms((sec_x - [obs_data.sec_x_obs])-mean(sec_x - [obs_data.sec_x_obs])), ...
    rms((sec_y - [obs_data.sec_y_obs])-mean(sec_A_y - [obs_data.sec_y_obs])), ...
    rms((sec_z- [obs_data.sec_z_obs]) - mean(sec_z- [obs_data.sec_z_obs])));
end;

%% print mean results


fprintf('A %5.2f %5.2f %5.2f\n', mean(sec_A_x - [obs_data.sec_x_obs]), ...
    mean(sec_A_y - [obs_data.sec_y_obs]), mean(sec_A_z- [obs_data.sec_z_obs]));

fprintf('B %5.2f %5.2f %5.2f\n', mean(sec_B_x - [obs_data.sec_x_obs]), ...
    mean(sec_B_y - [obs_data.sec_y_obs]), mean(sec_B_z- [obs_data.sec_z_obs]));

fprintf('C %5.2f %5.2f %5.2f\n', mean(sec_C_x - [obs_data.sec_x_obs]), ...
    mean(sec_C_y - [obs_data.sec_y_obs]), mean(sec_C_z- [obs_data.sec_z_obs]));

fprintf('D %5.2f %5.2f %5.2f\n', mean(sec_D_x - [obs_data.sec_x_obs]), ...
    mean(sec_D_y - [obs_data.sec_y_obs]), mean(sec_D_z- [obs_data.sec_z_obs]));

fprintf('E %5.2f %5.2f %5.2f\n', mean(sec_E_x - [obs_data.sec_x_obs]), ...
    mean(sec_E_y - [obs_data.sec_y_obs]), mean(sec_E_z- [obs_data.sec_z_obs]));

fprintf('F %5.2f %5.2f %5.2f\n', mean(sec_F_x - [obs_data.sec_x_obs]), ...
    mean(sec_F_y - [obs_data.sec_y_obs]), mean(sec_F_z- [obs_data.sec_z_obs]));

fprintf('GB %5.2f %5.2f %5.2f\n', mean(sec_GB_x - [obs_data.sec_x_obs]), ...
    mean(sec_GB_y - [obs_data.sec_y_obs]), mean(sec_GB_z- [obs_data.sec_z_obs]));

fprintf('HB %5.2f %5.2f %5.2f\n', mean(sec_HB_x - [obs_data.sec_x_obs]), ...
    mean(sec_HB_y - [obs_data.sec_y_obs]), mean(sec_HB_z- [obs_data.sec_z_obs]));

fprintf('I %5.2f %5.2f %5.2f\n', mean(sec_I_x - [obs_data.sec_x_obs]), ...
    mean(sec_I_y - [obs_data.sec_y_obs]), mean(sec_I_z- [obs_data.sec_z_obs]));

